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A SURVEY OF STRUCTURAL DYNAMICS 


OF SOLID PROPELLANT ROCKET MOTORS* 
by J. H. Baltrukonis 


ABSTRACT 


From the point of view of dynamics, a solid propellant rocket motor is an 
unique structure in that it is composed of a substantial mass of propellant 
material case-bonded to a relatively massless, thin-walled cylinder. The 
mechanical properties of the propellant are such that it contributes little to 
the stiffness of the composite^ structure but it does contribute to the dynamic 
characteristics of the structure because of its mass. Furthermore, due to the 
viscoelastic character of the propellant, it can be expected to provide con- 
siderable damping to the system. At this point in the development of the art 
of design of solid propellant rocket motors there are no clear-cut or well- 
founded methods to quantitatively evaluate the contributions of the propellant 
to the dynamic responses of the composite structure. It is clear that unless 
such methods are devised, it will be difficult to arrive at accurate missile 
designs . 

In the paper a survey of the dynamic problems of solid propellant rocket 
motors is presented starting from the simplest model thereof and proceeding, 
step-by-step, to the consideration of more sophisticated and realistic models. 
The consideration is restricted to infinitesimal deformations of propellant 
grains with linear mechanical properties. Substantial progress has been 
achieved towards the solution of many important dynamical problems. We shall 
attempt to summarize the pertinent developments and indicate along which lines 
we feel future study should proceed. A few illustrative solutions are included 
in areas wherein progress has been substantial. 

*An abbreviated form of this paper was presented at the International 
Conference on the Mechanics and Chemistry of Solid Propellants held at Purdue 
University, Lafayette, Indiana on April 19-21, 1965. 
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INTRODUCTION 


Structural dynamics concerns the analysis, by theoretical and/or experi- 
mental means, of the interactions of time-dependent loads and/or deformations 
externally applied to a structure or structural element and the internal 
stress and displacement response wherein inertial effects must be included in 
the analysis. It is the objective of this paper to present a survey of the 
field of structural dynamics of solid propellant rocket motors, to discuss 
those aspects of the subject which are of particular interest to the author 
and to recommend those areas in which further study should prove fruitful and 
rewarding. It is not our objective to present a bibliographical survey of the 
field and, consequently, many specific, individual contributions will probably 
be overlooked. There is no claim of uniqueness for this survey of the field 
nor do we maintain that it is complete since it is inevitable that a study of 
this sort will be biased to a large extent by the limitations, interests and 
viewpoint of the author. 

The first logical step in a task of this sort is to delimit the field 
under consideration. In general, we shall not be interested in the dynamics 
of complete rockets or missiles which may consist of several stages. Instead, 
we shall limit our concern to individual solid propellant rocket motors; i.e., 
casing plus propellant. We shall not be interested in any attachments to the 
casing such as rocket nozzles, guidance and control system assemblies, etc. 

It is not intended to imply that the dynamics of complete rockets is not im- 
portant, but rather the intention is to limit the scope of this presentation, 

A solid propellant rocket motor is an unusually complicated structure, 
at least from the standpoint of analysis thereof. It consists of a thin, 
circular cylindrical casing with domed end closure. The energetic propellant 
grain is bonded to the casing along its outer cylindrical surface. Frequent- 
ly, a rubber liner is interposed between the grain and casing for purposes of 
insulation. The flow of the gases developed by surface combustion of the solid 
propellant occurs through passages of relatively complex geometry within the 
grain to the one or more nozzles in the aft dome of the motor. The motor 
casing material is usually an elastic material such as steel though there has 
been a recent trend to casings wound of fiberglass filaments and impregnated 
with some sort of hard resin. 

The propellant material is a composite consisting of an elastomeric 
binder, a crystalline oxidizer and dispersed solid materials such as aluminum 
particles. This material is very compliant relative to the casing, is time- 
or frequency-dependent and is highly temperature-sensitive. Finally, the 
grain is very massive compared to the casing usually constituting from 80 to 
95 per cent of the total mass of the motor. 

Clearly, analysis of this structure for its responses to various dynamic 
stimuli is an imposing task indeed. The usual first step in the analysis of 
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• a complicated structure is the formulation of a tractable mathematical model. 
The structure involves several fundamental difficulties but probably the 
principal ones are (1) the complex geometry of the internal passages in the 
propellant grain; and (2) the fact that the motor is of finite length thereby 
introducing possible interaction of end effects. These two difficulties are 
of the same nature as the complications that have plagued elasticians from 
the very beginnings of the field of elasticity. Consequently, at the outset, 
it seems wise to idealize the structure to the extent that these complications 
are no longer present. Thus, we consider a mathematical model that is in- 
f initely-long with a circular perforation. We do not include a liner so that 
the motor consists only of grain and casing. We further restrict the present 
consideration to linear analysis; i.e., we consider only infinitesimal de- 
formations of linear materials. These restrictions may not be as severe as 
they seem when it is realized that most of the dynamic environments of interest 
will result in very small displacements within the linear range of the mater- 
ials under consideration. Materials properties investigations have revealed 
that, for small strains, propellants typically behave as linear ly-viscoelastic 
solids . 


ELASTIC GRAINS 


The mathematical model which we have thus far devised is still too com- 
plicated for initial studies. Thus, we introduce the further assumptions that 
the grain is perfectly elastic with no internal perforation and, since the 
casing is very stiff relative to the core, that the casing be perfectly rigid. 
Therefore, our initial, very crude model of a solid propellant rocket motor 
consists of an infinitely-long composite cylinder with a rigid outer layer 
and a very compliant, solid, elastic core. Investigation of such a model, 
however crude, will begin to yield valuable information concerning natural 
frequencies and displacement and stress fields corresponding to normal modes. 
Such information is of immediate utility, for example, in the design of the 
guidance systems whose reliability depends, to a large extent, on estimates 
of natural frequencies and mode shapes of the rocket motor. 

Other possible applications are dynamic loads analyses, staging studies, 
transportation and handling studies, etc. Over and above these practical 
applications, initial analyses of crude mathematical models provides the 
analyst with the experience required for subsequent refinement of the model. 
Additionally, a catalog of solutions is rapidly developed which may prove to 
be very useful as checks on the solutions of more refined models. 

The natural coordinate system for the formulation of the problem posed is 
the polar cylindrical system. There has long been activity within the general 
area of dynamics of elastic bodies in polar cylindrical coordinates. The 
earliest formulation of a problem by means of the dynamical equations of 
elasticity in cylindrical coordinates is due to Pochhammer (1) and Chree (2) 
who independently investigated longitudinal and transverse wave propagation 
in infinitely-long circular rods with traction-free surfaces. An excellent 
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survey of the history of this problem was presented by Abramson, Plass and 
Ripperger (3) but, nevertheless, it would be useful and interesting to mention 
a few of the more important investigations. Following Pochhammer and Chree 
the next study that should be mentioned is due to Ghosh (4) who derived dis- 
persion equations for longitudinal wave propagation in thick- and thin-walled, 
infinitely-long, hollow, circular cylinders with both surfaces free of trac- 
tion and with one surface traction-free and the other rigidly clamped. 

Bancroft (5) was among the first to publish numerical solutions of the 
dispersion equations for longitudinal wave propagation. In addition, he re- 
corded the dispersion equation for transverse wave propagation in infinitely- 
long, circular rods with traction-free surfaces but no numerical work was 
indicated. Hudson (6) carried out some numerical work with the dispersion 
equation for transverse wave propagation but arrived at the erroneous conclu- 
sion that there was only one mode of propagation of transverse waves. This 
error was subsequently propagated by Davies (7) and Kolsky (8), among others, 
and was finally pointed out by Abramson (9) who calculated several dispersion 
curves in the first mode of propagation of transverse waves. 

McFadden (10) was concerned with radial vibrations in hollow, thick- 
walled cylinders while Gazis (11) presented a study of plane strain vibrations 
in hollow cylinders with traction-free surfaces. The most complete study, to 
date, by means of the dynamical equations of elasticity, was presented by 
Gazis (12) on the longitudinal and transverse wave propagation and free vibra- 
tions in infinitely-long, thick-walled cylinders with traction-free surfaces. 
Other pertinent studies were reported by Herrmann and Mirsky (13), Mirsky and 
Herrmann (14, 15), Greenspon (16, 17, 18) Bird, Hart and McClure (19) and 
Bird (20). 

The field and constitutive equations for dynamic deformations of com- 
pressible elastic continua in polar cylindrical coordinates may be reduced 
to the following three equations of motion: 
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Solutions of these equations of motion are readily obtained by means of three 
displacement potentials as follows: 
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It may be verified by direct substitution that Eos. (1) are identically satis- 
fied by these components of displacement provided we take the displacement 
potentials as solutions of the following differential equations: 
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wherein c and c denote the dilatational and shear wave velocities, respec- 
tively, and are given by 


C 2 = k- 2 (g/p) 
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We recognize the differential equations (3) as wave equations, solutions of 
which are well known. On selection of appropriate solutions of these wave 
equations, the displacements follow from Eqs. (2) while the components of 
stress are obtained from the following stress-displacement relations: 
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Let us now illustrate the application of this general theory in the cal- 
culation of the natural frequencies and normal modes of the initial crude 
model for the solid propellant rocket motor. Since the grain is solid, the 
boundary conditions are given by 
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The simplest modes of vibration are the so-called axial-shear modes which 
are defined as those modes of free vibration in which the only non-zero 
component of displacement is that parallel to the axis of the cylinder. 
Furthermore, this axial displacement component depends only on the coordinates 
in the plane of a cross-section. These modes are of fundamental importance 
if for no other reason than that they are among the very few modes for which 
exact, closed-form solutions of the three-dimensional equations of elasticity 
are possible. There are, however, other reasons for concern with axial-shear 
modes. They represent limiting flexural or longitudinal modes; i.e., they are 
flexural or longitudinal wave modes with infinite wave-length. In accord with 
their definition we seek solutions of the equations of elasticity in the form 


u = u = 0 
r 6 


u z = W(r) e^ ut cos n0 


(7a) 

(7b) 


Substitution into Eqs. (1) results in the following single equaiton: 


W" + - W* 
r 
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wherein primes denote differentiation with respect to the argument. This is 
the well-known n^-order Bessel equation which has the following general 
solution: 

W(r) = C , J ( flr/a) + C _ Y„ (fl r/a) (9a) 

nl n n2 n 


wherein n is a dimensionless frequency coefficient defined by 


2 2 
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and C are arbitrary constants which must be evaluated such that the 
boundary conditions will be satisfied, and J and Y are the n^-order Bessel 
functions of the 1st and 2nd kinds, respectively. ?n view of Eqs. (7) the 
boundary conditions given by Eqs. (6) reduce to the following single condition 


W(a) = 0 



This condition will be identically satisfied provided that we take C to 
vanish and that 2n 


j (n) = o (io) 

n 

This result constitutes the frequency equation in the present problem. It 
defines a doubly-infinite family of natural frequencies and normal modes; 
i.e., corresponding to each value of n we find an infinite number of roots of 
Eq. (10). It is clear from Eq. (9a) that the n=0 modes are axisymmetric 
whereas the modes for which n=l are antisymmetric. In view of Eqs . (9a) and 
(7), we observe that all the stresses vanish with the exception of the shear 
stresses T and T and the latter vanishes for the axisymmetric modes. 
Clearly, tnis general type of vibration is pure shear in nature accounting, 
partially, for the name axialshear. These modes, among others, were thorough- 
ly investigated in Ref. (21). 

Another important class of vibrations occurs when the axial component of 
displacement vanishes identically yielding a transverse mode of vibration. 

Such a mode of vibration occurs when we take 

<J> = C . J (<s£) e* ,wt cos n0 (11a) 

nl n a 

ip = 0 (lib) 

x = c J (Si L ) e iwt sin n0 (11c) 

n2 n a 


It is immediately clear from Eqs. (2) that we obtain a plane strain mode of 
vibration wherein u z vanishes identically. Substitution into Eqs. (2) and, 
thence, into the boundary conditions given by Eqs. (6) results in the follow- 
ing frequency equation: 


J „-! <“> J n.l <' B > + J n + 1 <“> J „-l <*» > * 0 < 12) 

wherein we have used the recurrence relations for the Bessel functions. These 
transcendental frequency equations define a doubly infinite set of natural 
circular frequency coefficients which depend upon Poisson's ratio. In Refer- 
ence (22) this dependency is investigated in detail. The first four branches 
of the first four modes of vibration are plotted in function of < which de- 
pends only upon Poisson's ratio. Additionally, the displacement fields for 
several different modes are plotted. 
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The last problem of interest concerning this initial crude model is that 
of transverse wave propagation. In Ref. (23) this problem was treated in some 
detail. The dispersion equations were derived and several branches were 
plotted. Additionally, it was pointed out that the dispersion equations de- 
generate for infinite wavelengths into two uncoupled frequency equations de- 
fining the axial-shear and transverse vibrations modes which we have discussed 
above . 

Before proceeding to the refinement of our initial mathematical model, 
mention should be made of an interesting problem that arises in connection with 
Eq. (12), the frequency equation for free, transverse vibrations. It was men- 
tioned that the natural frequency depends upon Poisson's ratio. In Ref. (22) 
it was demonstrated that finite, real natural frequencies exist when the 
material of the core is ideally incompressible; i.e., when Poisson's ratio has 
the value 1/2. A question immediately arises: "How can natural frequencies 

exist for an incompressible material when it occupies the entire internal 
volume of the tank?" In Ref. (24) it was demonstrated that, as the core mater- 
ial tends to become incompressible; i.e., as Poisson's ratio tends to 1/2, 
the frequency spectrum tends to a simple line spectrum. Clearly, this kind of 
a frequency spectrum is physically impossible. That this should be the case 
is not surprising since the ideally incompressible material is a hypothetical 
material which cannot exist in nature. Nevertheless, the assumption of such 
a material is quite regularly used in practice. The behavior mentioned above 
is simply one difficulty that can arise from the application of such an as- 
sumption. Finally, we draw attention to a few interesting studies due to 
Magrab (25, 26, 27) concerning the displacement and stress fields in the solid 
grain due to steady-state, forced harmonic oscillation of the rigid case. 

On the basis of the studies cited we conclude that our initial crude 
model has been thoroughly investigated and that its steady-state response is 
adequately understood. Therefore, we proceed to a slight refinement of our 
mathematical model by allowing for a circular internal perforation of the 
grain. Now our model consists of an inf initely-long , two-layered cylinder 
with the outer layer ideally-rigid and the inner layer elastic. The boundary 
conditions in this case are given by the following: 



r=b r*b r=b 
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wherein a and b denote the external and internal radii, respectively, of the 
grain. Equations (13a) express the conditions that the internal grain perfor- 
ation be free of surface tractions while Eqs . (13b) result from the assumption 
of an ideal bond between the grain and rigid casing. 

For the axial-shear mode of vibrations we seek solutions in the form 
given by Eq. (7) so that Eq. (8) is the governing equation of motion while 
Eq. (9a) gives an admissible solution. In view of Eqs. (5), (7) and (9a) 
we find that all but two of the boundary conditions given by Eqs. (13) are 
trivially satisfied and these two conditions result in the following: 

(n — 0 | 1 | 2 , •• #) 

0 


C nl J n W * C n2 Y „ <«> * 0 


C „1 V (f! I > + C „2 V (a 7 > 


This is a homogeneous system of linear, algebraic equations in the unknown 
constants C nl and C n 2 * Such a system can have non-trivial solutions only if 
the determinant of the coefficients of the unknowns vanishes identically. 
Thus, we are led to the following frequency equations: 

Axisymmetric Modes: 

J (Q) Y. (ft £• ) - J. (0 - ) Y (ft) = 0 (14a) 

Q X cL x a o 

Anti-symmetric Modes: 


J (ft) Y * (ft T ) - J ’ (ft - ) Y (ft) = 0, n = 1, 2... (14b) 

n n a nan 

The zeros of these frequency equations were calculated in Ref. (21). 

In Ref. (28) the problems of plane strain vibrations and transverse wave 
propagation in the more refined model are treated. In the case of transverse 
wave propagation none of the boundary conditions given by Eqs. (13) is triv- 
ially satisfied so that a dispersion equation is obtained in the form of a 
6x6 determinant set equal to zero. The elements of the determinant are, in 
general, linear combinations of n t ^-order Bessel functions of the first and 
second kind. This dispersion equation is written in Ref. (28) and sample 
calculations were carried out in Ref. (29) for an incompressible grain. In 
Ref. (28) it is demonstrated that the dispersion equation governing transverse 
wave propagation degenerates for infinite wavelength to two uncoupled frequency 
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'equations. One of these is that for axial-shear modes of vibration as given 
by Eq. (14b) and the other is that governing plane strain vibrations. The 
latter frequency equation is a 4 x 4 determinant involving nth-order Bessel 
functions of the first and second kind set to zero. Four branches for each of 
the first four modes of vibration were calculated in Ref. (30) for an incom- 
pressible grain. ; 

The next refinement we can make in our mathematical model of the rocket 
motor is to allow for an elastic casing so that our new model consists in a 
composite cylinder of two concentric elastic layers both of which are infin- 
itely-long. In general, the stiffness of the grain will be considered small 
as compared with that of the casing material in deference to the origins of 
the problem . The understanding of the behavior of such a model should shed 
considerable light upon the rocket motor problem. To be sure, the outer layer 
which models the motor casing will usually be thin relative to its mean radius 
and, consequently, shell theory may be used in describing its behavior. How- 
ever, for the time-being, we shall treat both layers by means of elasticity 
theory as given in Eqs. (1) - (5). We shall return, subsequently, to the use 
of shell theory in treating the behavior of the outer layer. In order to as- 
sociate a given quantity with one or the other of the two layers, we append a 
1 when referring to the outer layer and a 2 when reference to the inner layer 
is intended. Thus, u Z 2 denotes the axial component of displacement in the 
inner layer while u z i denotes the same quantity in the outer layer. In accord 
with this convention, the boundary conditions are given by 
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For the axial-shear mode of vibration the only non-zero component of 
displacement is the axial component and we take it in the form 
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Direct substitution reveals that the equations of motion given by Eqs. (1) 
are identically satisfied for this displacement field. Substituting into 
Eqs. (5) we find that the only non-zero components of stress are the 
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following: 
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v T e see immediately that 8 of the twelve boundary conditions given by Eos. (15) 
are trivially satisfied and only the conditions given by Eqs. (15c,?,i,l) 
remain to be satisfied. Substitution into the latter conditions results in 
a homogeneous system of four linear, algebraic equations in the unknown 
constants C.-^ and C j 2 . The necessary and sufficient condition that non- 
trivial solutions of this system exist is that the determinant of the coef- 
ficients of the unknowns vanishes identic 
frequency equations: 
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wherein n = 0,1,2 ,.... On expansion and rearrangement, these frequency 
equations reduce to the follovnng : 
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wherein 
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The frequency equation given by Eq. (18) is an implicit function relating the 
density ratio P 2 /P 1 • the ratio of the shear modulii G 2 /G 1 » the radius ratios 
r^/a and r 2 /a and one or the other of the two frequency coefficients or 

^ 2 * i n Ref. (31) it is shown that this frequency equation includes Eqs. (10) 
and (14) as special cases, as it should. Additionally, various branches of 
the frequency equation are plotted in order to establish the conditions under 
which the earlier simpler solutions can be used. In general, it would seem 
that when the casing is reasonably thin, the assumption of a rigid casing is 
acceptable only for axisymmetric, axial-shear modes. 

For transverse (plane strain) vibrations of the composite cylinder we 
take the axial displacement to vanish identically and all displacements and 
stresses to be independent of the 2 coordinate variable. Under these condi- 
tions t rz and t z vanish identically and four of the twelve boundary condi- 
tions given by Eqs. (15) are trivially satisfied. Substitution into the re- 
maining boundary conditions results in a system of 8 homogeneous, linear, 
algebraic equations in the 8 unknown constants. The necessary and sufficient 
condition that non-trivial solutions of this system exist is that the deter- 
minant of the coefficients of the unknowns must vanish. Thus, the frequency 
equation is obtained in the form of an 8 x 8 determinant set to zero. In 
general, the elements of the determinant are linear combinations of n^- order 
Bessel functions of the first and second kinds. Calculations of natural fre- 
quencies have been carried out in Ref. (32) for a case wherein the outer layer 
(casing) is very much stiffer than the inner layer (grain). In this particu- 
lar case, the various branches of the frequency equation have been plotted and 
analyzed for modes wherein n = 1; i.e., for modes that have only a single nodal 
diameter. Modes of vibration have been identified that degenerate to pure 
grain and pure casing modes. A pure grain mode is defined as that mode of 
vibration existing in the grain when the case is perfectly rigid while a pure 
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casing mode occurs in an empty casing. It was shown that, for typical 
geometries and material parameter values, the lowest n = 1 casing mode is 
substantially higher than the lowest n = 1 grain mode. It is demonstrated 
that the natural frequency of the lowest casing mode is very sensitive to 
grain thickness. It may be possible to include this effect in an approximate 
solution based on shell theory wherein the mass of the grain is lumped with 
the mass of the thin casing and the stiffness of the grain is neglected. 

This possibility -should be exploited. Additionally, it was shown that the 
grain modes are relatively insensitive to variations in casing thickness 
indicating that, at least for the material parameters considered, the 
coupling of casing and grain rigidities is relatively weak. It follows that 
reasonable approximations to the natural frequencies of the grain modes can 
be obtained by assuming that the casing is perfectly rigid. The studies of 
Ref. (32) were limited to n = 1 modes. Work is in progress at The Catholic 
University of America to obtain similar results for n = 2,3,4 modes. The 
frequency equation has been programmed for machine computation, the program 
has been checked and data produced. It remains to analyze the data presently 
available and to supplement it as required. 

As previously mentioned, it is possible to use thin shell theory in 
treating the deformation of the casing rather than elasticity theory. However, 
since thin shell theory is not substantially simpler than elasticity theory in 
this particular application, there appears to be little profit especially 
since elasticity theory must be used for the grain. Sann and Shaffer (33) 
have performed an interesting study wherein shell theory was used for the 
casing. Although the frequency equations are developed for all modes of 
transverse vibrations, numerical results and analysis are presented only for 
axisymmetric modes. Two different solutions are presented: one applicable 

to a compressible grain and the other to an incompressible grain. It is 
found that with a compressible core, the axisymmetric mode has two uncoupled 
motions; one is a 'rigid' rotation of the casing with a twisting of the grain 
and the other is a purely radial motion. The latter motion is not present 
in an incompressible grain. Some simplified frequency equations are also 
presented for limiting extremes of rigidity and density ratios. 

The steady-state, transverse wave propagation problem for the two-layered, 
elastic cylinder remains for further treatment. None of the boundary 
conditions given by Eqs. (15) is trivially satisfied in this case so that the 
dispersion equation takes the form of a 12 x 12 determinant set equal to 
zero. This dispersion equation is readily written down but it remains to 
calculate its branches and perform the requisite analyses thereof. 

We are once again ready for an additional refinement of our mathematical 
model. There are two types of refinements that can be made, both of which 
are fraught with difficulties: (1) complicated internal perforation and 

(2) finite length. It would appear that each of these possible 
refinements requires individual treatment. 
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Let us first consider the problem of the complicated geometry of the in- 
ternal perforations of the grain. These internal passages are usually three- 
dimensional in character having shapes that are usually dictated by the con- 
siderations of internal ballistics. Little has been accomplished concerning 
the dynamics of such motors nor is it likely that a great deal will be achieved 
in the near future other than numerical solutions of specific problems. How- 
ever, it frequently happens that these internal grain perforations are two- 
dimensional in nature over substantial portions of the total length of the 
motor. In these portions the grain cross-section is circular, of course, with 
a star-shaped perforation. Therefore, it is not unreasonable to consider an 
infinitely- long grain with such a cross-section. A considerable volume of 
work has been accomplished with this model of a rocket motor and in the fol- 
lowing paragraphs we shall attempt to survey at least a portion of the work 
with which the author has a familiarity. 

For a mathematical model consisting of an infinitely-long, circular .grain 
with a complicated internal perforation, ideally bonded along its outer peri- 
phery to a rigid casing, the boundary conditions are that along the outer 
periphery the displacements must vanish while the periphery of the complicated 
internal perforation must be traction-free. The latter condition is expressed 
analytically as follows ; 


°r "r + T r9 n 9 = 0 


T r0 n r + Ve * 0 


x n + t. n. = 0 
rz r 0z 0 


along S(r,0) = 0 

( 20 ) 


wherein S(r, 0 ) = 0 defines the internal perforation and n r and n q denote 
the components in the radial and circumferential directions of the unit nor- 
mal drawn outwardly ^ith respect to the domain under consideration. In 
writing this condition we have taken n z , the axial component of the unit nor- 
mal, to vanish since the grain is cylindrical. Our initial considerations 
concerning this problem should be directed toward developing a method of solu- 
tion. For this reason we choose to simplify the problem by considering a 
solid cylindrical bar with its outer periphery having the same shape as the 
internal perforation of the grain. Thus, the boundary condition applicable 
to this bar problem is given exactly by £q. (20). Basically, the only dif- 
ference between the bar problem and the grain problem is that the displacement 
boundary conditions are absent in the former. Consequently, techniques suc- 
cessfully applied in the bar problem should also be successful for the grain 
problem. 
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Additionally, for the sake of simplicity, we concern ourselves only with 
the axial-shear mode of free vibrations. Thus, we seek solutions in the form 



<f> = X = o 

(21a) 


4> = V .Or, 8) e la>t 

(21b) 

Accordingly, Eqs. (2) reduce to 
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(22a) 


u = -e iut V, 2 T(r,6) 
2 1 w 

(22b) 

and Eqs. (3) degenerate to 

the following single Helmholtz equation: 



V ^ f = - P*— f 

1 G 

(23) 

Finally, as a consequence 
given by Eqs. (5) become 

of Eqs. (22), the stress-displacement equations 


°r = °e * “z = T re = 0 

(24a) 


2 iwt 

t = pur t — e 

rz K 3r 

(24b) 


2 1 3T iut 
T 0z = pa) F— 6 

(24c) 


Now, it can be shown* that the outward unit normal to S in the plane of S is 
given by 

-*■ f x -? VS 

n = ,n r t 

* See , for example , Wylie ( 34) 
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In view of Eqs. (24) and (25) the boundary conditions given by Eqs. (20) re- 
duce to the following single condition: 


3J13S ( 1 13S 

jr3r r 30 ; r 5T 


) 

. 



(26) 


Thus, the problem has been reduced to one of finding a solution of Eq. (23) 
such that the boundary condition given by Eq. (26) will be identically satis- 
fied. 


Closed form solution of the problem posed above is not feasible so we 
shall be content with approximate solutions. 

In Ref. (35) the collocation method is applied to solve the problem for 
a star-shaped boundary (with four star tips) given by 


S(r,6) = r - a - b cos 46 = 0 (27) 

We readily recognize the similarity of this family of plane curves to the 
boundary curve of the internal perforation of many common solid propellant 
rocket motors. Another advantage of this family is the fact that the circle 
is one curve of the family. Consequently, solutions can be degenerated to 
those for circular boundaries. Since the latter are available, we have a 
ready means for checking the results. The collocation method used consisted 
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t in taking a solution in the form of a finite sequence of solutions of Eq. (23) 
and satisfying the boundary condition given by Eq. (26) at a finite number of 
points on the boundary. Little is known concerning convergence of the col- 
location method applied to eigenvalue problems such as the present case. It 
is generally presumed that, provided a sufficient number of collocation points 
are used, the resulting eigenvalues will be reasonably accurate. Even less 
is known concerning the manner of distributing the collocation points along 
the boundary although it is generally presumed that the distribution becomes 
less important as the number of collocation points increases. In this study 
the first four natural frequency coefficients were calculated and plotted in 
function of the parameter b/a which governs the length of the star tips. 

When b/a = 0, the bar is circular and as b/a increases the star tip grows 
longer. The study concludes that, for the problem under consideration, the 
collocation method is very sensitive to the distribution of collocation points. 
Furthermore, little convergence is demonstrated for as many as seven colloca- 
tion points taken within an actant of the boundary. In view of the fact that 
it is not possible to obtain exact, closed form solutions in problems of this 
type, procedures for obtaining upper and lower bounds on the branches of the 
frequency equation are sorely needed. Only then can we be expected to make 
definitive statements concerning error in approximate procedures. Such 
bounding techniques frequently begin with estimates of the eigenvalues. Per- 
haps the value of the collocation method lies in its ability to provide these 
estimates fairly easily and quickly. 

Jain (36) has introduced a new kind of collocation procedure which elim- 
inates some of the difficulties of the boundary collocation method referred 
to above. He refers to the method as ’extremal point collocation’. He has 
applied the new method with striking success to the solution of several 
boundary value problems. The method requires the initial selection of a fin- 
ite number of collocation points. This selection is arbitrary as in boundary 
collocation. Instead of requiring that the error in satisfaction of the 
boundary conditions vanish at the collocation points , as in boundary collo- 
cation, in extremal point collocation it is required that the error at ad- 
jacent collocation points be equal in magnitude but opposite in sign. Fur- 
thermore, the error at the collocation points must be larger than that at any 
other boundary point. It is from the latter condition that the method de- 
rives its name. In order to satisfay these collocation conditions, an itera- 
tive procedure is required by means of which the distribution of collocation 
points is determined. Extremal point collocation has two distinct advan- 
tages. Selection of the collocation points is not arbitrary; instead it is 
determined by the method. In ordinary collocation there is no control over 
the maximum error or the distribution of error. In extremal point colloca- 
tion the error can nowhere exceed that at the collocation points and the 
error is uniformly distributed over the entire boundary. 

In Ref. (37) the problem defined by Eqs. (23), (26) and (27) was solved 
using extremal point collocation. In this study it is shown that the method 
is an effective technique for the calculation of eigenvalues. It is rela- 
tively simple to use provided that a large-scale computer is available. The 
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method .converges rapidly and yields reasonably accurate results even for a 
small number of collocation points. The question of accuracy requires fur- 
ther study by means of a method wherein the eigenvalues can be bounded. 

Let us return to the boundary conditions given by Eqs. (20) for further 
consideration. If we substitute into Eqs. (20) from Eqs. (24), we obtain 


( 


3r n r 


+ 


i at. 

r 39 n 6 


S 


0 


The quantity on the left is the normal derivative which defines the rate of 
change of ¥ in the direction of the normal to S. Thus, the boundary condi- 
tion that requires that the lateral surface of the bar be traction-free can 
also be written in the following form: 



S 


(28) 


This form suggests that there may be some advantage in formally mapping the 
bar cross-section onto a unit circle. If the conformal transformation is 
defined by 

w = w (5) 

wherein w = re defines the real plane while £ = Re 1 ' 1 defines the complex 
plane, the boundary condition given by Eq. (28) becomes 


i 




R=1 


(29) 


To be sure, satisfaction of this boundary condition is a trivial task com- 
pared to satisfaction of the boundary condition given by Eq. (28). However, 
we should hasten to point out that, by conformal transformation, we have 
simplified the task of satisfying the boundary conditions but we have com- 
plicated the task of finding solutions of the differential equation of motion 
since it, too, must be transformed. It can be shown* that the plane, 
Laplacian operator transforms according to 


3 2 T + 1 9T + 1 3 2 y 

9r 2 r 3r r^30 2 



(ill + 1 21 + 1_ 

3R 2 R 3R R 2 3>i 2 


( 30 ) 


*See, for example, p. 629 of Wylie (34) 
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and it follows, therefore, that, under the conformal transformation , Eq. (23) 
becomes 


3 2 V + 1 8? 1 3 2 <y _ pto 2 

3R 2 R 3r R 2 3ju 2 G 





(31) 


It is immediately clear that we will probably have difficulties finding solu- 
tions of this differential equation. Nevertheless, this formulation of the 
problem proves advantageous in certain circumstances. For example, common 
approximate methods of solution of problems of this type such as Rayleigh’s 
principle or the Ritz method require trial functions that satisfy the bound- 
ary conditions. In the transform plane such trial functions can be formulated 
with little difficulty. 


In Ref. (38) conformal transformation is used as outlined above in the 
solution of the problem of axial-shear vibrations of a star-shaped bar with 
a cross-section in the form of a four-lobed epitrochoid. The epitrochoidal 
boundary was chosen since its conformal mapping onto the unit circle is rela- 
tively simple and since it possesses the general characteristics of the solid 
propellant rocket grain perforation. To solve the problem in the complex 
plane the collocation method is applied wherein the boundary condition is 
identically satisfied and the error in satisfaction of the differential equa- 
tion is collocated at a finite number of points in the interior of the unit 
circle. The results were considered favorable though some difficulty was en- 
countered with regard to the spatial distribution of collocation points. 

Among the many methods available for the solution of eigenvalue problems 
the methods of Rayleigh and Ritz are probably the most familiar. These 
methods yield upper bounds on the eigenvalues but, in the absence of exact 
values, it is difficult to estimate the error in the approximate solution. 
Temple (39) suggested a method for estimating the error in each stage of an 
iteration procedure directed toward the precise determination of the lowest 
eigenvalue. Temple’s method can be interpreted as one which establishes the 
lower bound. A very powerful method for the determination of upper and lower 
bounds on all eigenvalues has been developed independently by Kohn (40) and 
Kato (41), as a generalization of Temple’s method. This method was applied 
in Ref. (42) for the calculation of the natural frequencies in axial-shear 
vibrations of circular and epitrochoidal bars. The resulting bounds were ex- 
tremely accurate for the fundamental natural frequency but the accuracy deter- 
iorated somewhat for the high natural frequencies. However, it should be 
pointed out that the bounds can be improved at will if one is willing to ex- 
pend the requisite additional effort. 

In an eigenvalue problem it is required to find one or more constants x, 
called eigenvalues, and corresponding functions g , called eigenfunctions, 
such that a differential equation 


M [8] = XN [6] 
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is satisfied throughout a domain D subject to certain boundary conditions on 
the boundary of D. In general, the domain D may be either a one- or a two- 
dimensional continuum. For many eigenvalue problems M [..] and N [..] are 
both linear, self-adjoint, positive-definite, defferential operators with the 
order of M greater than N. Under these conditions the eigenvalues are real 
and positive and the eigenfunctions form an orthogonal set. When the eigen- 
value does not appear in the boundary conditions, the eigenvalue problem is 
called special provided that the operator N has the form 


N [6] = gB 

where g is a prescribed continuous function which is positive throughout the 
domain D. Thus, the governing differential equation for special eigenvalue 
problems becomes 

M [B] = AgB 

We see immediately that Eqs. (23) and (31) have this form and, since the 
eigenvalue does not appear in the boundary condition, it is clear that the 
axial shear vibrations problem is a special eigenvalue problem in both the 
real and complex planes. Therefore, we can make use of the methods that have 
been developed for special eigenvalue problems. Collatz (43) has developed 
a procedure for obtaining upper and lower bounds in special eigenvalue prob- 
lems. Using a trial function which satisfies the boundary conditions, but 
not necessarily the differential equation, a few simple operations are per- 
formed and the upper and lower bounds result. However, these bounds may not 
be very good unless the trial function is close to the exact solution of the 
problem. The fact that the method is relatively simple to apply for any given 
function indicates that it might become much more useful if a systematic pro- 
cedure were added for ’’choosing" these trial functions. Such a procedure was 
developed by Appl and Zorowski (44). Another such method was developed in 
Ref. (45) and applied in the solution of the axial-shear vibrations problem 
for an epitrochoidal bar. The bounds obtained were adequate and subject to 
additional refinement but the effort involved is probably less than that re- 
quired to obtain Kohn-Kato bounds. 

We have discussed a number of techniques that should prove useful in the 
solution of the vibrations problem for the rocket motor model consisting of 
an infinitely- long, circular grain with a complicated internal perforation, 
ideally bonded along its outer periphery to a rigid casing. The more effec- 
tive of these techniques require trial functions that satisfy the boundary 
conditions. It is clear, therefore, that a key requirement is a technique 
for conformally transforming the circular domain with a complicated perfora- 
tion onto a circular annulus. Wilson (46) has developed such a technique 
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but it is limited to the case that the web fraction* is relatively small. 

For example, Arango (47) has shown that, when the perforation has four axes 
of symmetry, the error in mapping the external boundary increases rapidly 
when the web fraction increases beyond 0.5. When the web fraction reaches 
the value 0.9, the outer boundary has lost all semblance of a circle. This 
error arises due to the fact that Wilson treated the conformal transformation 
of an infinite domain with a hole into another such domain. Consequently, 
while the mapping function accurately transforms the internal boundary into 
the unit circle, the external boundary transforms only approximately into a 
circle. Rim and Stafford (48) have very recently presented a simple method 
of deriving approximate mapping functions in the form of low order polynomials 
which conformally transforms an annular region into one whose inner and outer 
boundaries are star-shaped and circular, respectively. The derivation is 
based on the Schwarz- Chris toffel transformation. This method has the same 
accuracy problems in transforming the outer boundary as Wilson’s method. 

This concern with the accuracy of the transformation of the external boundary 
is of substantial importance since, while the web fraction may be relatively 
small in the unburned solid propellant grain, it will increase toward unity 
as a consequence of the burning process. It should be clear, therefore, that 
it is essential to develop techniques for the conformal transformation of 
finite, doubly- connected domains. 

Laura (49) has developed such a technique based on numerical integration 
of a pair of dual integral equations derived by Kantorovich and Muratov (50) 
for the purpose of conformal transformation of an arbitrary , finite , doubly- 
connected region onto a circular annulus. Laura demonstrates that, if the 
domain under consideration has one or more axes of symmetry and one of the 
boundaries is a circle, the system of two integral equations simplifies con- 
siderably. Several illustrative transformations are presented including the 
transformation of a domain which is typical of a solid propellant rocket 
motor. The accuracy of the technique is exceptional. 

For purposes of illustration, let us consider a grain cross-section with 
4 axes of symmetry. An octant of this cross-section is shown in Fig. 1 be- 
tween the heavily-accented inner and outer boundaries. The outer boundary 
is a circle with a radius of 18 inches while the inner boundary has a maximum 
radius of 9-inches and a minimum radius of 3- inches. The inner and outer 
boundaries are further identified with the captions R = 1.0 and R = b = 2.61, 
respectively. A mapping function of the following form was used to con- 
formally transform (approximately but with more than adequate accuracy) the 
domain between the inner and outer boundaries in Fig. 1 onto a circular 


* Web fraction is defined as the ratio of the diameter of the circle circum- 
scribing the inner boundary of a doubly connected region to the diameter of 
the circle circumscribing the outer boundary. 
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Fig. 1 First approximation to a set of burning curves for a typical solid 
propellant rocket motor 
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(32) 


wherein p denotes the number of axes of symmetry. In this case we had 4 axes 

of symmetry and, therefore, we took p = 4. A twelve term mapping function 

was used; i.e., we took M = 1 and N = 10. A larger number of terms would 

have resulted in a more accurate mapping. The unknown coefficients, A • 

in Eq. (32) were evaluated by Laura’s method. Using this mapping function ^ 

seven additional circles in the complex plane, intermediate between the 

inner and outer boundaries, were transformed onto the curves shown in Fig. 1. 

Each of the contours shown is identified by the radius of the circle in the i 

complex plane upon which it maps. Thus, the curve R = 1.6 transforms into 
the circle with radius 1.6, the curve R = 2.0 transforms into the circle with 
radius 20, etc. If one were able to stop the burning process in a solid 
propellant rocket motor at seven instants of time between ignition and burn- 
out and to plot the shape of the inner boundary at each time, the so-called 
’burning curves’ would be obtained. The resulting contours would look much 
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like the curves of Fig. 1. Thus, as a first approximation, we regard the con- 
tours of Fig. 1 as a set of burning curves for a typical solid propellant 
rocket motor. Incidentally, it should be pointed out that, if actual burning 
curves were given, they could be conformally transformed onto circles in the 
complex plane with exceptional accuracy using Laura’s method. 

With the mapping function known it is not a difficult task to calculate 
the natural frequencies in the axial-shear mode for each of the regions of 
Fig. 1. The results would provide estimates of how the natural frequencies 
change during the burning process. There are any number of approximate pro- 
cedures available for this calculation. We choose to use Galerkin’s method. 
The applicable equation of motion is given by Eq. (31) which we choose to 
rewrite as follows : 


wherein 
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The quantity g(R,y ) was calculated using Eq. (32). The appropriate bound- 
ary conditions are 
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(34a) 


(34b) 


The first of these conditions expresses the fact that the displacement at the 
outer periphery of the grain must vanish since the grain is ideally bonded to 
a rigid case. The second condition above expresses the fact that the inner 
contour is free of tractions. The quantity a will take on 7 different values 
corresponding to the 7 different burning curves in Fig. 1. The following 
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trial function, satisfying the boundary conditions, was used: 
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We have chosen to ignore the y -dependency, which we are at liberty to do 
since we are attempting an approximation. If the y -dependency had been in- 
cluded, the resulting approximation would have been better than that obtained 
using Eq. (35). For these calculations we used a one-term trial function; 
i.e. , S = 1. Substitution from Eq. (35) into Eq. (33a) results in 

E (R,y) = f B n £w» (R) + ± (R) + X 2 g (R,y) W n (R)] (36) 

wherein e(R,y ) is an error function which does not, in general, vanish 
since, in general, W n is not an eigenfunction. Galerkin's method requires 
that the error function e (R, y) be orthogonal to the trial functions 
W n (R) throughout the domain of interest; i.e.. 


/ E(R, u) W (R) dD = 0, (n = 1,2, S) (37) 

D n 

Performing the indicated integration results in S homogeneous, linear, alge- 
graic equations in the S unknown constants B n . We obtain a frequency equa- 
tion by requiring that the determinant of the coefficients of the unknowns 
vanish identically. This frequency equation was solved 7 times for the lowest 
natural circular frequency coefficient corresponding to tlie 7 different values 
of a. The results have been plotted in Fig. 2. The horizontal coordinate in 
Fig. 2 is the web fraction but, since the web fraction will be a function of 
time as the burning process proceeds, the plot in Fig. 2 shows the variation 
of the lowest natural circular frequency coefficient with time throughout the 
burning process from ignition to burnout. We see that the natural frequency 
increases rapidly as the motor burns out which is not an unexpected result. 
Such information is of extreme importance in the design of a guidance loop for 
the vehicle in which the motor is to be used. 

It has been demonstrated that substantial progress has been made con- 
cerning the treatment of complicated grain perforations when these passages 
are substantially two-dimensional. Conformal transformation of the perfora- 
tion by Laura’s method should open the door to the use of a variety of accept- 
able methods for calculating natural frequencies in various modes of vibration 
and to the study of various steady-state wave propagation and forced vibra- 
tions problems. Incidentally, it is worth pointing out that Laura’s method 
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Fig. 2 Variation of the fundamental natural circular frequency coefficient 
in the axial-shear mode with web fraction for a typical solid propel- 
lant rocket motor 


is equally applicable to the study of problems wherein inertial effects are 
unimportant such as quasi-static internal pressurization. Thus far only 
axial-shear modes of vibration have been considered for grains with compli- 
cated perforations. Additional effort is required on transverse modes and in 
wave propagation studies. When the geometry of the perforation varies rapid- 
ly with distance parallel to the motor axis, application of these two-dimen- 
sional methods becomes questionable. Little, if any, work is being accom- 
plished on these three-dimensional problems. It would seem that future work 
in this area will be confined to the application of finite difference tech- 
niques. 

We come now to the consideration of the ultimate refinement of our model 
of a solid propellant rocket motor; consideration of a model of finite length. 
Most of the modern solid propellant rocket motors are short with length-to- 
diameter ratios of the order of unity not uncommon. Little has been accom- 
plished with the study of such dynamic problems. However, much is being done 
with quasi-static problems of this nature, particularly using finite differ- 
ence schemes of various types. The development of understanding of finite 
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length effects and the study of the validity of our two-dimensional results 
when applied to relatively short motors are probably the most important un- 
solved problems in solid propellant motor dynamics today. Some little work is 
being accomplished but there is a noticeable lack of literature on the subject. 
We should call attention to the experimental study of the vibrations of short , 
solid elastic cylinders due to McMahon (51). Extensive experiments were per- 
formed and comparisons with available theory presented. However, the avail- 
able theories were limited to thick disk and Timoshenko beam theories. It 
should be clear that much future effort should be expended in the area of 
finite length cylinders. 


VISCOELASTIC GRAINS 


The study of elastic grains has been reasonably active as can be judged 
from the detail of our discussion thereof. This is as it should be since, 
because of the elastic-viscoelastic correspondence principle,* solutions of 
elastic grain problems are intimately related to the solutions of viscoelastic 
grain problems. This principle states that the Laplace- or Fourier-transformed 
solution of a viscoelastic problem can be obtained from the Laplace- or 
Fourier- transformed solution of the associated elastic problem by replacing 
the elastic constants therein with appropriate functions of the Laplace or 
Fourier transform parameter. Therefore, in principle, it would seem that one 
can always obtain the viscoelastic solution associated with a given elastic 
solution. In practice, the situation is not quite so straight-forward since 
inversion of the viscoelastic solution remains to be accomplished. To be 
sure, exact inversion is always preferable when possible. In the practical 
situation this is rarely possible since the quantities replacing the elastic 
constants are usually measured material property functions of the Laplace or 
Fourier transform parameter available only in curve or tabular form. It 
should be clear, therefore, that a considerable effort should be mounted con- 
cerning approximate techniques for performing Laplace and Fourier transforma- 
tions and inversions. In this regard we should mention the methods discussed 
by Schapery (56) and Cost (57) for the Laplace transform and by Solodovnikov 
(58) and Papoulis (59) for the Fourier transform. Assuming that approximate 
inversion methods will be available, our principal concern should be with 
elastic problems since these could be readily converted to viscoelastic solu- 
tions. Clearly, therefore, those difficulties that arise in elastic problems 
will remain with us in viscoelastic problems; for example, the finite length 
difficulty is still troublesome. Summarizing, it seems evident that visco- 
elastic problems present no fundamental difficulties that are not inherent in 
the associated elastic problems; the computational problem is more complicated 
and tedious but the basic qualities of the two types of problems remain the 
same. 


*See Alfrey (52), Lee (53) and Sips (54, 55 ) 
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Because of the intimate relationship of elastic and viscoelastic problems 
and because of our rather lengthy discussion of elastic grain problems, we 
shall not belabor the viscoelastic problem area. However, we will mention a 
few pertinent investigations and present an illustrative solution, just enough 
to impart the flavor of the problem. 

Gottenberg (60) has reported the results of an experimental investigation 
of steady-state, transverse vibrations of a long, steel cylindrical tube con- 
taining an inert propellant with a circular perforation. Many bending modes 
were detected and the resonant frequencies and mode shapes compared with the 
predictions of a Timoshenko beam theory in which the bending stiffness of the 
propellant was neglected relative to the casing stiffness but the additional 
mass of the inert propellant was included. The comparisons were quite adequate 
for engineering purposes. Substantial difficulties were encountered in detect- 
ing modes other than bending modes. However, one axisymmetric , longitudinal 
mode and a few breathing modes were identified. No theory was available for 
comparison with these modes. 

Henry and Freudenthal (61) reported the results of an extensive analytical 
investigation of the steady-state, forced vibrations of a thick-walled visco- 
elastic cylinder contained by and bonded to a thin, cylindrical shell. Only 
axisymmetric solutions were considered. Complex frequency response functions 
were determined which may be easily used for arbitrary and random inputs by 
means of the well-established methods of generalized harmonic analysis. The 
study is broken down into three steps: (1) the thick-walled cylinder, (2) the 

thin shell, and (3) the composite cylinder. For a thick-walled, elastic 
cylinder the solutions of Eqs. (1) - (5) are developed in the usual manner ex- 
panding the normal tractions on the inner boundary and the normal and tangen- 
tial tractions on the outer boundary in Fourier trigonometric series in the 
axial coordinate. The boundary conditions on the lateral surfaces are identi- 
cally satisfied but, rather than requiring the axial normal and sheer stresses 
to vanish over the ends, the axial normal stress and the radial displacement 
are required to vanish. As a consequence of this relaxation of the boundary 
conditions, the physical significance of the results has been obscured to a 
certain extent. The difference between this ’finite-length 1 solution and the 
infinite-length solution is not clear. In order to obtain the viscoelastic 
solution, the elastic-viscoelastic correspondence principle is invoked wherein, 
for a forced vibrations problem, the elastic constants are replaced by complex 
material properties functions of frequency. Membrane theory is used for the 
thin elastic casing and a higher-order shell theory is included in an appendix. 
The latter should be used when the membrane theory is inadequate; for example, 
for external loads of acoustic nature. The assumption that the propellant and 
case form a continuous structure at their interface requires continuity of dis- 
placements and, therefore, produces strong coupling of the motions of grain and 
casing. This coupling is considered of primary importance and, therefore, the 
interaction is treated rigorously. The internal pressure is taken as harmonic 
and a large volume of numerical results is presented for the fundamental radial 
mode. It is concluded "that any analysis of a solid fuel rocket motor which 
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does not take into consideration the viscoelastic effects of the propellant 
would not give a true overall picture". Another important observation is ar- 
rived at by varying the grain geometry in simulation of the state of the com- 
posite structure at various stages of the burning process. The results indi- 
cate that, for certain ranges of frequency, there occurs a considerable in- 
crease in the amplitudes of stress and displacements. The implication is that 
the critical period in rocket operation would occur in the later stages of the 
burning process when the propellant is almost completely burned out. Thus, 
the important frequencies might be close to the fundamental radial mode of the 
casing . 

In Ref. (62) the stress response to pressure transients has been investi- 
gated in an infinitely-long, two-layered cylinder having a thin, incompressible 
elastic outer layer and an inner layer of an incompressible, two-parameter 
Voigt material. This composite structure was taken as a crude mathematical 
model of a solid propellant rocket motor. The particular problem of interest 
concerned the circumferential stress response of the case to the pressure 
transients induced at ignition in the grain perforation. Two different pres- 
sure programs were studied in some detail: a square ignition pulse and a 

triangular ignition pulse followed by a pressure, constant in time. By solving 
the elastic problem first and then invoking the correspondence principle to 
obtain the viscoelastic solution, it was shown that the stress response is 
rather insensitive to the shape of the pulse. However, it was also demonstrated 
that the duration of the transient is important. When the duration of the 
transient is an order of magnitude smaller than the natural period of the cyl- 
inder in the radial mode, the stress response barely reflects the presence of 
the ignition pulse. However, as the duration of the transient approaches the 
natural period, the circumferential stress approaches the stress level corres- 
ponding to the pressure magnitude of the transient. The implication as far as 
casing and igniter design should be evident. If one is to economize on casing 
weight but nonetheless maintain rapid and positive ignition, the igniter must 
be designed such that the pressure transient will persist for a time no larger 
than one order of magnitude smaller than the natural period of the motor in the 
radial mode. 

Lockett (63) presented the results of a study of significant importance with 
regard to the analysis of transient responses in viscoelastic materials. He 
discussed the effect produced by a rapid, but not discontinuous, change in pres- 
sure at the surface of a spherical cavity in an infinite, viscoelastic medium. 
Invoking the correspondence principle, Lockett obtains the solution to the 
viscoelastic problem from the Fourier transform of the associated elastic 
problem by replacing therein the elastic shear and bulk moduli with the visco- 
elastic complex shear and bulk modulus functions of circular frequency. Thus, 
he is subsequently concerned with performing inverse transformations. It is 
the manner of performing these inversions that is significant. Before proceed- 
ing to the actual inversion procedure, Lockett discusses the nature of modulus 
functions that should be used and the character of the pressure- time history. 
Most solutions to viscoelastic problems appearing in the literature are either 
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. left in the integral form or are specialized to some simple spring-dashpot 
representation for the complex moduli. However, it has been shown by Kolsky 
and Shi ( 64 ) that, in general, these simple models do not adequately represent 
the behavior of real, viscoelastic materials. Thus, if only the solution to 
a particular problem is required, it would seem advisable to use the experi- 
mental data directly in the numerical inversion of the Fourier-transformed 
solutions. If it is desired to keep a number of parameters in the solution, 
then mathematical models should be chosen which fit the experimental data well, 
even though they may not correspond to a simple spring-dashpot configuration. 

The simple forms of modulus functions corresponding to spring-dashpot models 
only have an advantage in the solution of simple problems when it may be pos- 
sible to evaluate the integrals explicitly. In his subsequent numerical work 
Lockett takes a constant bulk modulus; i.e., he assums that the material is 
elastic in dilatation, and for the complex shear modulus he uses an analytical 
expression that has been derived from curve fitting to experimental measurements 
performed on real, viscoelastic materials. For a pressure-time history Lockett 
selects a time function which has the following complex Fourier transform: 


0, 0 < to 

P (u>) = ^ P q / iu), _< a: _< u) 2 

0, u 2 — w — < ’° 


( 38 ) 


wherein the range ( 05) defines the range in which accurate measurements 

are available of the complex shear modulus. It follows that, since the Fourier 
transform of the pressure-time history is a factor in the transformed solution, 
the exact solution in both the elastic and viscoelastic problems is a finite 
integral over the range ( 05) within which the shear modulus is accurately 

known. The following pressure-time history is obtained by inversion of the 
Fourier transform given by Eq. ( 38 ): 
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where Si(x) denotes the sine integral which has been tabulated by Jahnke and 
Emde ( 65 ) and Abramowitz and Stegun (66), among others. For “ 2 very much 
larger than w^this pressure program has the shape shown in Fig. 3 . Lockett 
considers this loading in its own right and not as an approximation to an 
exactly square loading. The latter is not obtainable in practice, so the load- 
ing shown in Fig. 3 is more realistic. Obviously, the rise-time of the pres- 
sure program can be made shorter by increasing 00 2 up to the upper limit of the 
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Fig. 3 First version of the alternate pressure program 


frequency range in which the shear modulus is well-defined. The Fourier inver- 
sions of the transformed solutions are now performed numerically since integra- 
tion is over a finite frequency range only. Lockett completes the solution and 
arrives at some interesting conclusions concerning the differences between 
elastic and viscoelastic solutions. 


Making use of Lockett's suggestions the problem of dynamic internal pres- 
surization was solved in Ref. (67) for an infinitely-long , incompressible, 
viscoelastic cylinder case-bonded to a thin elastic tank. The pressure-time 
history of Fig. 3 has some of the character of a pressure-time program of a 
solid propellant rocket motor but it has two important deficiencies: the ap- 

plied pressure is negative for time less than zero and the pressure rise starts 
at negative time. Both of these deficiencies were corrected by appropriate 

shifting of axes so that Eq. (39) becomes 


P(t) 

P 

o 



to. 


Si (to^t — it ) — Si ( w^t - tr~ 


(40) 


This expression defines the one -parameter , family of curves shown in Fig. 

4. We see that these curves possess most of the characteristics of actual 
pressure-time programs for solid rocket motors. We see a finite rise time 
which is controlled principally by the value of to 2 and given very closely by 


t R = 27t/ “2 


(41) 
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Fig. 4 Pressure-time program applied to the inner cylindrical surface of the 
core 


The family of curves rises to a maximum and then. oscillates about a quasi- 
static curve whose shape depends upon the value of the ratio nij/ In 

Ref. (67) a value of 1000 was used for this ratio and for this value the os- 
cillation is about a curve that is approximately parallel to the dimensionless 
time axis at least for a substantial period of time. Since attention was 
focused on the immediate dynamic effects of the pressurization it mattered 
little that the pressure program decreases gradually for long times. The 
Fourier transform of the pressure program is given by 
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This transform is non-zero over a limited frequency range only so that, ulti- 
mately, the inversion integrals can also be evaluated by numerical methods. 

The shear storage and loss moduli used in the investigation were the actu- 
ally measured data shown in Figs. 5. These material properties were introduced 
into the transformed solution by expressing the shear storage modulus in the 
following form: 
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wherein 6,3,62,62, > %> Tj_ j Tg j , t m and N are arbitrary constants 

which must be evaluated such that Eq. (43a) will^yield a good approximation to 
the shear storage modulus shown in Fig. 5(a) in the frequency range of interest. 
Equation (43a) is the analytical representation for the shear storage modulus 
for the (2N+1) -parameter Maxwell material. Thus, use of this expression implies 
that we have idealized the solid propellant as a linearly-viscoelastic, (2N+1)- 
parameter Maxwell material. It is clear in Eq. (43a) that G Q constitutes the 
shear storage modulus at zero frequency and, on the basis of extrapolation of 
the data of Fig. 3(a), the value selected was 

6 q = 400 psi 

The remaining parameters in Eq. (43a) were evaluated by means of a method in- 
troduced by Schapery (68) wherein the T n 's were arbitrarily assigned and the 
C n 's evaluated by solving a system of N linear algebraic equations obtained by 
the collocation method. Five collocation points were used and the resulting 
analytical expression agreed with the curve of Fig. 5(a) to within the width 
of a pencil line. The analytical expression for the shear loss modulus for 
the (2N+1) -parameter Maxwell material, which is also required, is given by 

N WT 

G" (») = l G 2-s- (43b) 

n=l n 1 + 


It is clear, therefore, that the shear storage and loss moduli are not indepen- 
dent since the unknown parameters were evaluated from shear storage modulus 
data and these same parameters allow immediate calculation of the shear loss 
modulus by means of Eq. (43b). This result is not surprising and was recog- 
nized by Gross (69) when he calculated the following exact inter-relation be- 
tween these two material property functions: 

G" (to) =| J“ 6' (a) 2 “ -g- da (44) 

a - a) 


As an interesting experiment, the parameters calculated by means of Schapery' s 
method using Eq. (43a) and the shear storage modulus data of Fig. 5(a) were 
used in Eq. (43b) to calculate the shear loss modulus. The results are plotted 
as a dotted curve in Fig. 5(b). We observe a very substantial disagreement 
between this result and the measured experimental data. This discrepancy has 
not been explained. Rather than using the experimental data of Fig. 5(b), 

Eq. (43b) was used for the shear loss modulus. 
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Fig. 5(a) Dynamic shear storage modulus as a function of frequency for 
Hercules Powder Co. CYH solid propellant 
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Fig. 5(b) Dynamic shear loss tangent as a function of frequency for Hercules 
Powder Co. CYH solid propellant 

Making use of the transformed pressure program given by Eq. (42) and the 
material property functions shown in Eqs. (43) the transformed stresses and 
displacements in the two-layered cylinder were calculated. The inversion in- 
tegrals were evaluated numerically with relative ease since they were inte- 
grals over a finite frequency range only. The results for the circumferential 
stress at the internal perforation are shown in Figs. 6. These results display 
some interesting features to which we should draw attention. First, we draw 
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attention to the small perturbations displayed by all solutions for time less 
than zero and for time greater than zero although, in the latter time regime, 
they may be obscured by larger oscillations. These perturbations of the 
calculated responses are clearly due to the small perturbations present in 
the pressurization program. These small oscillations can be ignored since 
our interest concerns the gross responses due to the gross pressure increase. 
Next, quasi-static solutions corresponding to all transient solutions were 
obtained (and plotted in Figs. 6) by taking the specific weights of both case 
and grain equal to zero. Now, we observe that, for the larger perforation 
(7.0-incher) the circumferential stress at the internal perforation is 
compressive despite the fact that the accompanying circumferential strains 
are tensile. See Figs. 6(a). This same effect was noted in Reference (70) 
and was explained therein. We summarize the explanation herein for purposes 
of completeness. The circumferential stress may be thought of as consisting 
of two components: a hydrostatic compressive component due to the compression 
of the core against the restraint offered by the case and a tensile 
component due to radial growth of the core under pressurization. If the 
core is sufficiently thin, such that the case restraint is important, 
then the compressive component is larger than the tensile component and the 
resulting circumferential stress is compressive even though the circumferential 
strain is tensile. This is undoubtedly the situation when a = 7.0 inches. 

It is clear that, with the appropriate combination of material properties 
and geometry, the opposite could also be true. As a matter of fact, we 
see from Fig. 6(b) that, for a = 1.3 inches and for both rise times, the 
circumferential stress at the internal perforation and the accompanying 
circumferential strain are both tensile. It seems obvious in this situation 
that case restraint plays the subordinate role. 

Now, we observe that for the 40 millisecond rise time the transient and 
quasi-static solutions coincided, for all practical purposes, for all 
times and in all cases considered. It seems clear that the 40 millisecond 
rise time is long compared to the fundamental natural period of the 
cylinder and, under these conditions, the transient and quasi-static 
formulations of the problem are equally valid thus accounting for the good 
agreement. However, for the shorter rise time we observe a substantial 
transient at small times and for the larger perforation. See Figs. 6(a). 

We also observe that the transient does not occur for the shorter rise 
time and for the smaller perforation. There appear to be two possible 
explanations of this behavior. For the larger perforation we have shown 
that the elastic casing plays a major role. The presence of this elastic 
storage element could account for initiation of the transient oscillation 
and for sustaining it thereafter. For the smaller perforation we have 
already demonstrated that the case plays a subordinate role and since there 
is a lack of the elastic storage element, the transient oscillation, if 
initiated at all, would attenuate very rapidly as a consequence of the 
dissipative quality of the core material which plays the major role. There 
is another candidate explanation. In Fig. 6(a) for the smaller rise time 
we see that the period of the transient oscillation, which is, of course, 
the fundamental natural period of the composite cylinder, has a value of 
about 4 or 5 milliseconds. Thus, since the rise time very nearly coincides 
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with the natural period of the cylinder, we expect to excite a substantial 
oscillation which will rapidly attenuate. For the smaller perforation, the 
natural period changes and the rise time no longer coincides and we do "not 
expect to see a large transient oscillation. It is difficult, without 
extensive study, to choose between these explanations although we do favor the 
latter. There are, undoubtedly other features of these results that bear 
further disscussion but we feel that we have drawn attention to the salient 
features . 

Achenbach (71, 72) has recently completed two investigations of general 
import. The first of these concerns the transient, dynamic response of an 
incompressible, elastic grain with an ablating circular port case-bonded to a 
thin, elastic case. Despite the fact that this study concerns an elastic 
grain, it is discussed here because of its close relationship with the 
second study. An internal pressure program in the form of a step function in 
time is applied and an approximate solution obtained by asymptotic integration 
of the equation of motion. It is demonstrated that burning time affects the 
period of the vibratory response but not the amplitude thereof. In the 
second investigation. Ref. (72), the dynamic response is calculated in an 
infinitely-long, composite cylinder consisting of a thich-walled cylinder of 
an incompressible, viscoelastic material case-bonded to a thin elastic shell. 
Two dynamic loadings are treated: a time step in internal pressure and a time 
step in external pressure. The latter loading condition is applied to an 
inoperative motor in a missile launching silo. The problem is solved, 
numerial results plotted and discussed for a grain of a standard, linear 
solid (three-parameter Maxwell material.) The manner of solution is discussed 
when measured data is used for the grain shear modulus. 

We can judge by the character of these last few references that a great 
deal has been accomplished concerning transient, dynamic responses in 
viscoelastic grains. However, additional effort is required to perfect and 
simplify the application of the approximate methods for performing Fourier 
transforms and inversions. Then we can look forward to solving many transient 
dynamic problems using as input measured pressure programs and material 
property functions. 

A great deal has been achieved in the solution of general wave propagation 
problems in linear viscoelasticity. See, for example. References (8), (73) - 
(91). This list of references is by no means complete nor is it a listing of 
the more important papers. Despite the fact that so much effort has been 
expended in investigating this general problem, thereby testifying to its 
importance, I was unable to locate a single reference involving a solution of 
direct applicability to solid propellant rocket motors. It seems evident 
therefore, that a considerable effort is required to elucidate those 
environments that give rise to wave propagation problems and in solving the 
problems posed. 
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CONCLUSION 


We have completed a survey of dynamic problems in solid propellant 
rocket motors. Attention was restricted to linear analysis of elastic 
and viscoelastic grains. It was not intended that this survey should 
be a bibliographical study but rather a presentation of progress in the 
field in an attempt to establish the current state-of-the-art. Many 
important related areas have been omitted in the interests of conservation 
of time and space. Consider, for example, nonlinear dynamic problems, 
shock problems, coupled thermo-mechanical problems, acoustic instability 
problems, and experimental measurement of propellant properties. Effort 
should be continued and expanded in these important areas among others. 

We have accomplished a great deal with regard to the dynamic problems of 
solid propellant rocket motors, however, additional effort is absolutely 
essential to extend the goals that we have already reached. 
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NOTATION 


r,0,z 


radial, circumferential and axial coordinate variables 
of polar, cylindrical coordinate system 


t 


time 


P 


mass density 


G 


shear modulus 


u ,u„,u 


A 

V 2 


<M,x 

C 

C 

c 

s 

K 


radial, circumferential and axial 
displacement 

du u 1 

cubical dilatation = — — 

or r r 

Laplacian differential operator = 

ai_ 

displacement potential functions 

dilatational wave velocity = (G/ptc 2 ) -1 "'^ 

1/2 

shear wave velocity = (G/p) 

ratio of shear wave velocity to the dilatational 
wave velocity 

radial, circumferential and axial components of 
normal stress 


components of 


du 

+ z 


du t 

xr 

di + ii + i a 2 + 

dr 2 r dr r 2 d0 2 


1 /9 


x * x .t shear stress components 

r6* 6z* rz 

X Lame f s second elastic constant = 2Gv/(l-2v); 

dimensionless circular frequency coefficient 


a,b t r ,r 2 radii 

w, Q circular frequency and circular frequency coefficient, 

resp. 


n 


mode number 



Bessel functions of the first and second kinds, resp., 
and of the n-th order 


hi 


i 


> 


(-D 


1/2 


Primes over a variable denote ordinary derivatives of the variable with 
respect to its argument; e.g. * f f (r) = (d/dr)f(r) 
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